% plot GPD
% Need simu.mat, variable_wage_XX.mat, nump_XX.mat
clearvars;

load simu.mat nump dlogP;
nump_base=nump;
dlogP_base=dlogP;
clear nump dlogP;
load variable_wage_01.mat nump dlogP dlogP_sub;
dlogP_g01=dlogP;
dlogP_subg01=dlogP_sub;
clear dlogP dlogP_sub;
load nump_01.mat nump
nump_g01=nump;
clear nump;
load variable_wage_05.mat nump dlogP dlogP_sub;
dlogP_g05=dlogP;
dlogP_subg05=dlogP_sub;
clear dlogP dlogP_sub;
load nump_05.mat nump
nump_g05=nump;
clear nump;
load variable_wage_00.mat nump dlogP dlogP_sub;
dlogP_g00=dlogP;
dlogP_subg00=dlogP_sub;
clear dlogP dlogP_sub;
load nump_00.mat nump
nump_g00=nump;
clear nump;

dlogP_ng01 = [dlogP_g01; dlogP_subg01(:,2)];
dlogP_ng05 = [dlogP_g05; dlogP_subg05(:,2)];
dlogP_ng00 = [dlogP_g00];%; dlogP_subg00(:,2)];

figure(10)
loglog(sort(dlogP_base(dlogP_base>0))/40000,(length(dlogP_base(dlogP_base>0)):-1:1)/length(dlogP_base(dlogP_base>0)),'-',...
sort(dlogP_ng00(dlogP_ng00>0))/40000,(length(dlogP_ng00(dlogP_ng00>0)):-1:1)/length(dlogP_ng00(dlogP_ng00>0)),':',...
sort(dlogP_ng01(dlogP_ng01>0))/40000,(length(dlogP_ng01(dlogP_ng01>0)):-1:1)/length(dlogP_ng01(dlogP_ng01>0)),'--',...
sort(dlogP_ng05(dlogP_ng05>0))/40000,(length(dlogP_ng05(dlogP_ng05>0)):-1:1)/length(dlogP_ng05(dlogP_ng05>0)), '-.',...
'LineWidth',2);
legend('Baseline','\gamma_m=0','\gamma_m=0.1','\gamma_m=0.5');
xlabel('Avalanche size dlogP conditional on dlogP>0');
ylabel('Relative frequency');
ylim([10^(-7) 10^0]);
grid; grid minor;
set(gca,'FontSize',20);

figure(11)
loglog(sort(nump_base(nump_base>1))/40000,(length(nump_base(nump_base>1)):-1:1)/length(nump_base(nump_base>1)), '-',...
sort(nump_g00(nump_g00>1))/40000,(length(nump_g00(nump_g00>1)):-1:1)/length(nump_g00(nump_g00>1)), ':',...
sort(nump_g01(nump_g01>1))/40000,(length(nump_g01(nump_g01>1)):-1:1)/length(nump_g01(nump_g01>1)), '--',...
sort(nump_g05(nump_g05>1))/40000,(length(nump_g05(nump_g05>1)):-1:1)/length(nump_g05(nump_g05>1)) ,'-.',...
'LineWidth',2);
legend('Baseline','\gamma_m=0','\gamma_m=0.1','\gamma_m=0.5');
xlabel('Avalanche size L/N');
ylabel('Relative frequency');
xlim([10^(-4.5) 10^0]);
ylim([10^(-7) 10^0]);
grid; grid minor;
set(gca,'FontSize',20);

disp('frequency of L=\pm 1')
disp([(sum(nump_base(:,2)==1)+sum(nump_base(:,2)==-1))/length(nump_base(:,2)) ...
(sum(nump_g00==1)+sum(nump_g00==-1))/length(nump_g00) ...
(sum(nump_g01==1)+sum(nump_g01==-1))/length(nump_g01) ...
(sum(nump_g05==1)+sum(nump_g05==-1))/length(nump_g05)]);

disp('V(L)/E(L)');
Lpos_base=nump_base(nump_base(:,2)>1, 2)-1; % subtracting the Calvo
Lpos_g00=nump_g00(nump_g00(:)>1)-1; % subtracting the Calvo
Lpos_g01=nump_g01(nump_g01(:)>1)-1; % subtracting the Calvo
Lpos_g05=nump_g05(nump_g05(:)>1)-1; % subtracting the Calvo
disp([var(Lpos_base)/mean(Lpos_base) ...
var(Lpos_g00)/mean(Lpos_g00) ...
var(Lpos_g01)/mean(Lpos_g01) ...
var(Lpos_g05)/mean(Lpos_g05)]);